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Abstract. We present a numerical and analytical study of the thermal fragmen- 

o 

tation of a turbulent flow of interstellar hydrogen. We first present the different 
O i 1 dynamical processes and the large range of spatial (and temporal) scales that need 

to be adequately represented in numerical simulations. Next, we present bidimen- 
' sionnal simulations of turbulent converging flows which induce the dynamical con- 

densation of the warm neutral phase into the cold phase. We then analyse the cold 
structures and the fraction of unstable gas in each simulation, paying particular 
attention to the influence of the degree of turbulence. When the flow is very tur- 

u : 

bulent a large fraction of the gas remains in the thermally unstable domain. This 
unstable gas forms a filamentary network. We show that the fraction of thermally 
unstable gas is strongly correlated with the level of turbulence of the flow. We 
then develop a semi-analytical model to explain the origin of this unstable gas. 
This simple model is able to reproduce quantitatively the fraction of unstable gas 
observed in the simulations and its correlation with turbulence. Finally, we stress 
the fact that even when the flow is very turbulent and in spite of the fact that 
a large fraction of the gas is maintained dynamically in the thermally unstable 
domain, the classical picture of a 2-phase medium with stiff thermal fronts and 
local pressure equilibrium turns out to be still relevant in the vicinity of the cold 
structures. 
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1. Introduction 

The dynamics of the interstellar gas is of great importance to understand the formation of 
structures such as molecular clouds and it is therefore determinant in the star formation 
process. 

In this respect, the neutral atomic hydrogen (HI) is crucial since it represents more 
than 50% of the ISM and since the molecular clouds form through the condensation of 
HI gas. Because of the numerical stiffness of the problem, previous numerical models 
attempting to simulate the ISM at a scale of about lkpc and to form molecular clouds 
self-consistently have not considered heating and cooling functions that allow thermal 
bistability between 100 and 8000 K (c.g Vazquez-Scmadeni et al. 1995, 1996, Korpi et 
al. 1999, Ballesteros-Paredes et al. 1999) or have a numerical resolution which is not 
appropriate to adequately describe this physics down to the scale of the CNM structures 
(Gazol et al. 2001). Therefore it is currently unclear and indeed almost unexplored to 
what extent the physics of HI may or may not have a significant influence on the formation 
and the evolution of molecular clouds. 

From the theoretical point of view (Field et al. 1969, McKee & Ostriker 1977, Wolfirc 
et al. 1995) as well as from the observational one (Low et al. 1984, Boulanger & Perault 
1988, Kulkarni & Heiles 1987, Troland & Hcilcs 2003), it is now well established that HI 
is a thermally bistable medium which at thermal equilibrium and for a thermal pressure 
close to about (in the vicinity of the sun) 4000 K cm -3 , can be in two different thermo- 
dynamical states, namely a warm and diffuse phase (WNM) and a cold and dense one 
(CNM) roughly in pressure equilibrium. 

The linear stability analysis (Field 1965), the quasi-static thermal front propagation 
(Zel'dovich & Pikel'ner 1969, Penston & Brown 1970) and more generally the non-linear 
development of a single structure (see Meerson 1996 for a review and Sanchez-Salcedo 
et al. 2002 for a recent study) have been under investigation for a long time and are 
reasonably well understood. However, it is only recently that the behaviour of a thermally 
bistable flow in the fully dynamical or turbulent regime has been investigated. 

1.1. Previous work 

In the context of solar chromosphere and corona, Dahlburg et al. (1987) and Karpen et 
al. (1988) perform 2D simulations of a gas that may undergo thermal instability. They 
have considered random velocity perturbations in an initially thermally unstable gas 
and study the development of cold structures. They note that large amplitude velocity 
perturbations can prevent the formation of condensed structures. 

Send offprint requests to: E. Audit, P. Hennebelle 
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Kovalenko & Shchekinov (1999) and Hcnncbcllc & Pcrault (1999) have simultane- 
ously investigated the possibility that a converging flow may dynamically trigger the 
thermal transition from the WNM (thermally stable) phase into the CNM phase. They 
perform ID simulations and show that if the perturbation lasts long enough (more than 
a cooling time) and is strong enough (velocity must be comparable to the sound speed), 
the thermal transition is possible, i.e part of the WNM phase condenses in CNM leaving 
a cold structures embedded in a warm surrounding medium. Their underlying idea is 
that an external forcing like bubble expansion or any phenomena generating systematic 
or turbulent motions may promote the formation of cold structures. A similar picture has 
been investigated by Koyama & Inutsuka (2000) who simulate a shock propagating in 
HI and include H 2 formation. Hennebelle & Pcrault (2000) have also considered the case 
of a magnetic field, important in the context of the ISM, which is initially oblique with 
respect to the velocity field. They show that whatever the value of the magnetic field 
the thermal condensation is still possible provided the initial angle between B and V is 
small enough. The smallest angle below which thermal condensation is always possible, 
is about 20°. Due to magnetic tension the condensation is more difficult for intermediate 
magnetic intensity than for a stronger field. 

This paradigm has been further investigated by Koyama & Inutsuka (2002) who 
simulate in 2D a shock propagating in HI. They show that several structures of CNM 
form close to the shock interface and find that the velocity dispersion of the cold structures 
is about 5-6 km/s, i.e a fraction of the sound speed of the WNM. 

Gazol et al. (2001) have performed 2D simulations of HI at a scale of about 1 kpc and 
a numerical resolution of about 5 pc with a turbulent forcing that mimics star formation. 
They found an important fraction (~ 50%) of thermally unstable gas in their simulation. 

Kritsuk & Norman (2002a, 2002b) have performed 3D simulations of HI with a ther- 
mal forcing that mimics the random fluctuations of the heating rate derived by Parravano 
et al. (2002). Their aim is to show that interstellar turbulence may be generated by ther- 
mal instability. 

Finally, Piontek & Ostriker (2004) have recently performed 2D simulations with the 
aim of studying the development of the magneto-rotational instability as well as the 
thermal instability in the magnetised warm atomic interstellar gas. 

1.2. Outline of the paper 

In this paper we further investigate the dynamical properties of HI gas by the mean of 
2D numerical simulations. We explore the behaviour of HI focusing on the formation 
of cold structures during the collision between two streams (converging flow) of WNM. 
We put special attention to the structure morphology, their internal velocity dispersion 
and to the fraction of gas in the cold and warm phase and in an intermediate thermally 
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unstable state. Our aim is to understand how those features depend on the turbulence 
of the flow. 

In the second section of the paper we present the equations, describe the thermal 
processes, the numerical scheme and the initial and boundary conditions. We also discuss 
the drastic numerical resolution which is needed in order to properly describe the HI flow. 
In the third section we present our numerical results for a large scale flow which is weakly 
or very turbulent and present statistical analysis of the simulations. In the fourth section 
we emphasize the correlation between turbulence and the fraction of thermally unstable 
gas. We then develop a semi-analytical model to understand the physical mechanisms 
which are responsible for this phenomena. The fifth section summarises our results and 
concludes the paper. 

2. Equations, numerical requirement and initial conditions 

2.1. Equations and notations 

In this paper we consider the Euler equations for an optically thin gas. The gas is able to 
cool radiatively and is heated up by an external radiation field. The equations governing 
the evolution of the fluid are the classical equations of hydrodynamics, where a cooling 
function is added in the energy conservation equation: 

d t p +V.H -0 (1) 

d t pu + V.[pu®u + P]=Q (2) 
d t E + V.[u{E + P)] =-C(p,T) (3) 

where p is the mass density u the velocity, P the pressure, E the total energy and C the 
cooling function (see next section). The gas is assumed to be a perfect gas with 7 = 5/3 
and with a mean molecular weight p — I Amu, where tuh is the mass of the proton. 

The above equations are solved using a second-order Godunov method for the con- 
servative part. The cooling is applied after the hydrodynamical step using an implicit 
scheme and subcycling when the cooling time is much smaller than the time step. 

2.2. Thermal processes 

In order to simulate HI it is important to adopt a model for thermal processes which on 
one hand is realistic enough and on the other hand not too computationally expensive. 
We have followed closely the work of Wolfire et al. (1995, 2003), trying to include only 
the processes which are the most relevant to our study. The UV flux is taken equal to 
Go/ 1.7, where Go is the so called Draine's flux (Draine 1978). 

We have only kept the following cooling processes, which are dominant in the physical 
conditions of our simulations : 
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- cooling by fine-structure lines of CII (92 K) 

- cooling by fine-structure lines of 01 (228 and 326 K) 

- cooling by H (Lya line) 

- cooling by electron recombination onto positively charged grains 

The heating process is the photo-electric effect on small grains and polyaromatic 
hydrocarbons due to the far-ultraviolet galactic radiation. We have not taken the heating 
due to cosmic rays and to the soft X-rays into account since these contributions appear to 
be small. Finally in order to calculate the ionisation we use the approximation proposed 
by Wolfire et al. (2003). 

We have compared each term with the results given in Wolfire et al. (95) as well as 
the thermal equilibrium curve. Very similar results have been obtained. 

2.3. Resolution requirement 

The thermal condensation of HI gas involves 4 different spatial scales, each of them 
related to a different physical mechanism that we describe below. 

2.3.1. First scale: cooling length in the WNM 

The cooling length of the WNM, A coo i is defined as the product of the cooling time, 
Tcooi = C V T / '(Cp), by the sound speed, C s . It represents the scale at which the WNM 
is non-lincarly unstable (Hennebelle & Perault 1999) , i.e velocity fluctuations of spatial 
extension ~ A coo i and amplitude ~ C' s , can undergo a thermal contraction from WNM 
phase into CNM phase. For typical WNM parameters, n wnm ~ 0.5 cm" 3 , T wnm — 8000 
K ? A coo i and is about 10-20 pc whereas for a higher density of n wnm ~ 3 cm" 3 , reached 
for example in shock or dynamically compressed layer A coo i is smaller and about 1-3 pc. 

2.3.2. Second scale: fragment size 

The second important scale of the problem is the typical size of the fragments of CNM. 
It is given by the cooling length divided by the density ratio between the CNM and the 
WNM, ( ~ 100. This is due to the fact that a piece of initially warm gas contracts until 
it reaches thermal equilibrium. Therefore A strU ct = A coo i/C — 0.1 pc. We would like to 
point out the fact that resolving this scale is essential in order to properly describe the 
structure of the cold gas. Note that A str uct is the typical size of the structures in the 
direction along which the gas has condensed. In the other directions, if the collapse is 
anisotropic, the scale is a priori different and is rather related to the initial scale of the 
WNM fluctuation which has condensed into the CNM structure. 
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2.3.3. Third scale: the Field's length 

The third scale is the Field length (Field 1965) or conduction length, which is also the 
typical size of the front between the 2 phases. It is given by Apicid = y/ (k(T)T)/C 
In the WNM, the Field's length is about 1CT 1 pc and is about 1(T 3 pc in the CNM. 
In order to have an accurate description of the thermal front it is necessary to include 
non-ideal processes, such as viscosity and thermal diffusion. Since the front propagation 
velocity is proportional to AFicid (Zel'dovitch & Pikcl'ner 1969, Penston & Brown 1970), 
the consequences of not resolving properly this scale would be first to overestimate the 
propagation of the thermal front and second to overestimate the fraction of thermally 
unstable gas located in the front. We thus conclude that in order to simulate properly the 
HI gas it is essential that the effective Field length (close to the numerical resolution) in 
the simulation is small compared to A str uct ■ This ensures that the growth of structures due 
to heat diffusion remains small compared to the growth of structures due to dynamical 
processes and that the fraction of thermally unstable gas stabilised by heat diffusion and 
located in thermal fronts is small compared to the fraction of cold gas. The importance 
of resolving the Field length has been recently carefully analysed by Koyama & Inutsuka 
(2004). 

2.3.4. Fourth scale: shocked CNM 

The fourth scale is a consequence of dynamical processes. The internal sound speed of the 
CNM is about 10 times smaller than the sound speed of the WNM. Consequently, it is 
expected that HI structures experiment supersonic motions with a typical Mach number, 
M, equal to about 10. This may occur through the collision between 2 fragments of CNM 
or when a fragment of CNM is swept up by a shock wave. Assuming that the polytropic 
index of the CNM is about 1 (i.e. CNM is isothermal), this leads to a density increase of 
Pshock = M 2 p cnm . The typical size of the shocked layer is thus: X s hock = ^struct / M 2 which 
is about ~ 10~ 3 pc for M = 10. If the numerical resolution is larger than X s hock, the 
largest density reached in the simulation is underestimated and the size of the shocked 
layer is artificially broadened. 

As we have seen, simulating the thermal condensation in HI satisfactorily, would 
require to treat spatial scales ranging from few 10 pc to about ~ 10~ 3 pc simultaneously. 

2.4. Boundaries and initial conditions 

Due to the large range of spatial scales involved in the interstellar medium it is not 
possible to properly resolve the scales associated with the microphysics (see previous 
section) and at the same time to treat the large scale flow of the WNM. It is therefore 
necessary to compromise between the different scales. 
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Fig. 1. The figure illustrates the inflowing velocity field for e = 0.5 (dotted line), e = 2 
(dashed line), e = 4 (solid line) and e = 6 (dot-dashed line). As can be seen, we have 
kept the same phase for the 3 curves. 

Since we are interested in the formation of small dense structures, we have chosen to 
resolve the small scales and to use a simple prescription to model the large scale velocity 
field of the WNM. Therefore, the simulations used in this paper are 2D simulations on a 
1000 2 grid. In order to solve the cooling length the size of the box is L = 20pc leading to 
a numerical resolution of 0.02 pc. This numerical resolution ensures that \ s truct is well 
resolved, and that it is larger than the effective Field's length (equal to about one pixel). 
However with this resolution a shocked structure is not well described for Mach numbers 
larger than about 3. This means that the largest densities reached during supersonic 
collisions are underestimated in our simulations. 

The size of the simulation box, 20 pc, is marginal in the sense that it is not enough 
to describe the evolution of a large scale structure of WNM. The large scale flow is 
then imposed by choosing boundary conditions that mimic a converging flow at large 
scale. The upper and lower sides of the simulations have free boundary conditions while 
on the left (resp. right) side the gas is injected with a velocity Vi n (y) (resp. —Vi n {y)). 
The density and the pressure of the inflowing gas are chosen such that the gas is at 
thermal equilibrium in the branch of the WNM phase. It is thus thermally stable when 
it penetrates in the simulation box. 

The function V in describes the velocity of the inflowing gas and is given by: 

v in (y) = v e -((^/ 2 )/o.6Lr {1 + u{y)) (4 ) 

where U(y) = e ^ Cj cos(kjy + a,) 

This velocity field represents a stream of gas with a transverse spatial extension of about 
L/2 and an average velocity equal to Vq. In order to study the influence of turbulence, 
this field can be modulated by the function U(y). e is the amplitude of this modulation; 
fcj = 2ir/\i is the wave number and the wave-length A, lies between 2 cells and L/4. We 
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Fig. 2. Density and velocity fields in the central region (4x4 pc around the center 
located at x — 10 pc y = 10 pc) of the simulation box (20x20 pc 2 ) for e = 0.5 (weakly 
turbulent flow) at time t=18.93 Myr. The longest arrow represents a velocity of about 
15 km/s. 



have chosen a power law spectrum of index —1 for the modulation. Therefore the Ci are 
defined by Ci oc fer 1 and ^2 c i = 1- The on are random phases which lie between and 
2n. If e is small then the flow remains essentially laminar whereas it becomes turbulent if 
e is greater than 1. Fig. ^ illustrates the shape of the inflowing velocity field for different 
values of e (e=0.5, 2 and 4). 

The gas injected into the simulation box is at thermal equilibrium in the WNM branch 
and has a constant pressure equal to P = 710~ 13 erg cm~ 3 . The gas temperature and 
density are respectively ~7100 K and n = 0.76 cm -3 . 
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Fig. 3. Gas mass fraction in the pressure-density diagram for the e = 0.5 simula- 
tion. The colours correspond respectively to gas mass fraction (in arbitrary units) of 
(1,5,10,50,100,200,1000) for (yellow, orange, red, green, dark blue, light blue, pink). The 
full black line is the thermal equilibrium curve. The green lines are the isothermal curve 
T=5000 K and T=200 K and the blue line is the Hugoniot curve of shocked gas corre- 
sponding to our initial conditions. 

We have compared ID results involving a thermal transition between the 2 phases 
with the code used by Hennebelle & Perault (1999) and obtained very similar results 
within few percents of accuracy even though the code used in this work does not include 
viscosity and thermal conduction. 

3. Results 

In this section we qualitatively describe results for a weakly turbulent (e = 0.5) and for 
a very turbulent forcing (e = 4). We then present a statistical analysis for the 4 cases, 
e = 0.5, 2, 4 and 6 considering the gas fraction in the different phases and the properties 
of the CNM structures formed in the simulations. 

Initially, there is only warm gas in the simulation box. Then the two facing incoming 
flows generate a region of higher density and pressure in the center. This region is ther- 
mally unstable and cold structures start to form. After some time (from 5 to 15 Myr), the 
simulation reaches a permanent regime. That is the mass fraction in the different phases, 
the statistical properties of the cold structures, their abundance etc... remain constant. 
All the results presented below are given for this permanent regime. 

3.1. Weakly turbulent flow 

When the flow is weakly turbulent, the ram pressure of the 2 incoming flows induces the 
formation of a high pressure central layer which is bounded by 2 accretion shocks. This 
high pressure gas has a density of about 3 cm -3 and a temperature of about 10 4 K. Due 
to the high pressure in the central region a roughly homologous expanding velocity field 
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density and velocity fields, t= 18.93 My 
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Fig. 4. Density and velocity fields for a very turbulent forcing (e — 4). The longest arrow 
represents a velocity of about 17 km/s. 

develops in the transverse direction and part of the gas escapes the box continuously. The 
velocity field as well as the density field in the central region can be seen in Fig. [3 which 
displays a zoom of 4x4 pc 2 around the box centre (located at x — 10 pc, y = 10 pc). 
The layer is out of thermal equilibrium and its central part becomes thermally unstable. 
Part of the warm gas condenses into cold gas. This behaviour is indeed very similar to 
the ID situation studied by Hennebelle & Perault (1999). However in the present case 
because of the weak turbulence, the cold gas layer is unstable and fragments in several 
parts as can be seen in Fig. |3 where about 10 fragments have formed. Their size ranges 
from about 0.05 to 0.3 pc. Very sharp thermal fronts (2-3 pixels) bound the structures 
and connect them to the warm surrounding medium. Due to the high pressure induced 
by the incoming flow, the density of the structures is rather high and ranges from a few 
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100 to ~ 2000 cm -3 . The fragments are sometimes connected by a thin layer of low 
density gas. 

This situation is very similar to the 2D numerical simulations of Koyama & Inutsuka 
(2002) who study the thermal transition induced by a shock propagating in the warm 
phase. As in our case, Koyama & Inutsuka (2002) find that the shocked layer fragments 
in few cold clouds of about 0.1 pc. 

Fig. |3 gives a more accurate description of the thermal state of the gas in the simula- 
tion. It shows the gas fraction in the pressure-density diagram. The full line is the thermal 
equilibrium curve. The green lines are the isothermal curves T=5000 K and T=200 K 
and the blue line is the Hugoniot curve of shocked gas corresponding to our initial con- 
ditions. Most of the warm gas is located between the Hugoniot curve and the isothermal 
curve T=5000 K whereas most of cold gas is very close to the thermal equilibrium curve. 
There is almost no thermally unstable gas. 

3.2. Very turbulent flow 

Fig. 0] displays the density and velocity fields in the whole box for the very turbulent 
forcing (e = 4) at time t=18.93 Myr. The stiff shear of the gas that penetrates in the 
box generates turbulence quickly. 

The interface between the 2 flows is very irregular and significantly distorted. The 
density field is much more complex than in the previous case. As illustrated in Fig. [7J 
less material is at a density of more than few 100 cm~ 3 and a significant fraction of 
the gas is at intermediate density around 5 cm -3 and temperature smaller than 5000 
K, i.e in a thermally unstable state. Moreover, as can be seen by comparing Figs. [3] and 
the average thermal pressure is higher in the weakly turbulent case. This is due to 
the fact that the kinetic energy of the incoming flow (which is almost constant is all 
simulations) is converted into both thermal pressure and turbulent motions. Therefore, 
if the turbulent motions are high, the thermal pressure should be lower. 

As can be seen in Fig. [S] and which displays respectively the density and velocity 
fields and the temperature field and Mach number of Fig. 01 between x = 8 and 12 pc 
and y = 10.5 and 14.5 pc, the thermally unstable gas is very filamentary and presents a 
complex structure (We use the word filament for the elongated structures seen in our 2D 
simulations. In 3D these could become either sheets or filaments, though we believe that 
filaments are more likely since sheets would be more easily broken by turbulent motions) . 

Several interesting features appear. The different phases are highly interwoven with 
pockets of warm gas embedded into filaments of cooler gas. This is particularly obvious 
in Fig. H3 around x ~11.5 pc and j/~11.5 pc. 
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density and velocity fields, t= 18.9" 
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Fig. 5. Spatial zoom of Fig.0] The longest arrow represents a velocity of about 13 km/s. 



In spite of the presence of this unstable gas, the sharp fronts separating the cold 
and warm phases are still obvious as can be seen for the structure located at x =9 pc 
and y =13.8 pc. This relatively unsurprising result means that even in a very dynamical 
medium the 2-phase behaviour is not suppressed and may indeed be locally rather similar 
to the standard equilibrium 2-phase model. 

Sharp transitions can be seen between the warm phase and the filaments of thermally 
unstable gas as well (see for example the front at x ~ 11.2 , y ~ 11.3 pc). 

There are more cold structures than in the previous case (e = 0.5) but they are rela- 
tively less dense. This suggests that turbulence promotes the fragmentation of thermally 
unstable medium and that due to the more random motions (i.e less organised than for 
the case e = 0.5) the average thermal pressure (due to the high ram pressure) of the 
medium is reduced. 
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Fig. 6. Temperature and Mach number corresponding to Fig. |21 The longest arrow rep- 
resents a mach number of about 10. 
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Fig. 7. Same as Fig. for e = 4. 
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The cold structures are clearly associated with the unstable gas since they are often 
linked to other cold structures by a filament of unstable gas and sometimes embedded in 
such a filament. 

Note that the size of the smallest cold structures is close to the mesh of our 
simulation. This means, as pointed out in sect. [2 that numerical convergence is clearly 
not reached for those small structures. Moreover as can be seen in Fig. |SJ the cold 
fragments have a complex internal structure which is also smoothed out by the lack of 
resolution. 

We would like to stress the fact that in the dynamical regime, the 2-phase behaviour 
is clearly not erased and that the description of the medium as a continuum of phases 
would not be correct. There is thermally unstable gas but it is highly structured, very 
filamentary and connected to the denser thermally stable gas which is very differently 
structured. This means that new phenomena due to the thermal nature of the flow rather 
than simple disappearance of the two phase behaviour occur in this regime. 

3.3. Gas fraction in the different phases 

In this section, we analyse the temperature, density and pressure distributions obtained 
in the 2 cases presented previously (e = 0.5,4) and 2 additional cases (e = 2 and 6). 
For simplicity and in order to give simple trends, we call thermally unstable gas the gas 
having temperature between 5000 and 200 K, whereas the gas with a temperature below 
200 K is defined as cold gas. 

Fig. [S] shows histograms of the fraction of gas (X) as a function of temperature, 
density and pressure in the four cases e=0.5, 2, 4 and 6 for 4 time steps: t = 11.6, 15.8, 
17.9 and 20 Myr. 

For the case e = 0.5, the distributions at the 4 time steps presented are very similar, 
which means that the permanent regime is reached quickly. The fraction of thermally 
unstable gas is low (~ 10%) and the fraction of CNM is about 30%. There is almost 
no gas at intermediate density (n ~ 10 cm -3 ) whereas X is almost constant for density 
between 100 and 1000 cm -3 . It then drops stiffly for n > 1000 cm -3 . The pressure 
fluctuates over 1 order of magnitude. The peak at log(P/P me d) — —0.1 is due to the 
preshocked gas whereas the peak at ~ 0.5 is the postshocked gas (P mc d is the median 
pressure) . 

In the case e = 2, the fraction of cold gas is slightly lower (20%) than for e = 0.5 
whereas the fraction of thermally unstable gas is larger (20%). The fraction of gas at 
density above 10 cm~ 3 is lower at the first time step than later on which means that the 
permanent regime is longer to reach. 
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Fig. 8. Fraction of function of temperature, density and pressure for e = 0.5, e = 

2, e = 4 and e = 6 at times t=11.6, 15.6 and 20 Myr. 



For the cases e = 4 and 6, these trends are even clearer. It is seen that the distribution 
at time t = 11.6 Myr is significantly different from the distributions at time t = 15.8, 
17.9 and 20 Myr. The fraction of dense gas is smaller at t = 11.6 Myr by approximately 
a factor 3. This means that the permanent regime is much longer to reach than for the 
previous cases. Anticipating the mechanism that will be discussed in Sect. 01 we interpret 
this result as the consequence of the fact that turbulence is able to stabilise the thermally 
unstable gas making it able to last longer and consequently delaying the formation of 
the cold gas. In the same way, it is seen that the fraction of cold gas when statistical 
equilibrium is reached (<=17.9 and 20 Myr) is about 10% and thus significantly lower 
than for the cases e = 0.5 and 2. On the contrary, the fraction of thermally unstable gas 
is higher (about 30%) and the drop at intermediate density (~ 10 cm~ 3 ) is less clear. 
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3.4. Analysis of the CNM structures 
3.4.1. The analysis 

We now investigate the properties of the CNM structures in more details. We compute 
their surface and aspect ratio, their internal and relative (to each other) velocity dis- 
persion, as well as internal density, temperature and pressure distributions. A precise 
definition of a " structure" is needed. We define a structure as a connex set of pixels hav- 
ing a density above a given threshold, n s = 30 cm~ 3 . This definition has a clear physical 
meaning since the front that connects the 2 phases is stiff and therefore the structures 
do not depend very much on the arbitrary threshold n s . 

Once having identified the structures, we compute the inertia matrix / defined by 
I xx = J px 2 dxdy, I yy = J py 2 dxdy and I xy — I yx — J pxydxdy. It admits the 2 eigen- 
values Ai > A2. The aspect ratio r is then defined by r = ^/-WAi. We also consider 
the velocity of the structure, V s t rU ct, its surface, S, average density, < p >, temperature, 
< T >, and internal dispersion velocity defined as: 



where N is the number of pixels in the structure. 

Last, we compute the temperature, density and pressure variance for a structure 
defined as: 



3.4.2. Morphology and properties of the CNM 

Fig. displays the distribution of the density, temperature, pressure and velocity of the 
structures as well as their variance (see Sect. 13. 4T1 for the definitions). It also shows their 
surfaces (in pixels) and their aspect ratio. In order to obtain the distributions, we have 
added the structures obtained in 20 different time steps (separated by about 0.5 Myr in 
time) leading to a total number of about 800 to 1500 structures of CNM (depending on 
the value of e). 

The density histograms confirm the trend mentioned in the previous section that the 
structure are much denser in the weakly turbulent case than in the turbulent one by 
approximately a factor 10. They are consequently colder (factor 3) and have a higher 
thermal pressure (factor 3) for e = 0.5 than for e = 4. 

The density and temperature variances (respectively ~0.6 and ~0.7 for e = 4 and 
~0.7 and ~1 for e = 0.5) indicate that these quantities vary significantly inside one 
structure around their mean value. The pressure variance is significantly lower (~ 0.3) 
which indicates that the structures are on average not far from mechanical equilibrium. 




(5) 




(6) 
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Fig. 9. Properties of CNM structures for the 3 cases e = 0.5 (dashed line), e = 2 (dot- 
dashed line) and e = 4 (full line). The histograms display the density, temperature, 
pressure, velocity, distributions as well as their variance, the size (in pixels) and the 
aspect ratio distribution. See text for definitions of these quantities. 



The structure velocity ranges from to 8 km/s and is slightly lower in the weakly 
turbulent case (e = 0.5) than in the turbulent one. This is consistent with the results 
obtained by Koyoma & Inutsuka (2002) who found a comparable velocity dispersion 
for their clouds. The internal velocity dispersion is also slightly lower in the weakly 
turbulent case, for which it peaks at ~ 0.35 km/s, whereas it peaks at ~ 0.45 km/s for 
e = 4. This indicates that most of the internal motions are subsonic with respect to the 
internal sound speed (~0.8 km/s) of the CNM structures and is consistent with the small 
pressure variance. 
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Fig. 10. Evolution of the semi-analytical model for a model with (9, £) = (0, 0) (bottom 
panel) and for the same model with (9, S) = (0, 0.8/r coo j) (Upper panel). The dashed 
lines are the 5000 K and 200 K isothermal curves. 

Finally, it is seen that the mean surface of the structure is about 20-30 pixels wich 
gives a typical length of about ~ 5 pixels or ~0.1 pc. As already mentioned the structures 
are smaller when turbulence is higher. The structures have a mean aspect ratio of about 
0.5 for e = 4 and 0.35 for e = 0.5. 

4. Turbulence and thermally unstable gas 

In this section we present a more detailed analysis of the role of turbulence in producing 
the filaments of thermally unstable gas discussed in the previous section. We first present 
a semi-analytical model that describes the transition of a gas clump from the warm to 
the cold phase and propose an explanation for the origin of the unstable gas. We then 
compare the predictions of the model with the numerical results. 

4.1. A semi-analytical model for the evolution of a fluid particle 
4.1.1. Formalism and physical interpretations 

In order to understand how turbulence generates thermally unstable gas, we have devel- 
oped a semi-analytical model which is presented in the appendix. The underlying physical 
idea of the model is simply to follow a fluid element and to calculate the evolution of its 
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density and temperature by computing the cooling and heating of the gas and its geo- 
metrical evolution. It is based on a Lagrangian form of Eqs. Q-Q and the assumptions 
that the fluid element can be simply described by its semi-minor and major axis length 
a(t) and b(t), its internal pressure P(t), the external pressure P cx t and its density p(t). 
With these assumptions the fluid element can be described by the following equations: 

P+ pO =0, (7) 

£+ P/ P 6 =--C{p,T), (8) 
P 

In these equations 9 is the divergence of the velocity field, £ = 0.5(d x u x — d y u y ) 
(where the axis x and y are the eigenframe of the shear tensor, see appendix), a(i) = 
ao exp(fy(9/2 + Ti)dt) and b(t) = aoexp(J*(9/2 — T,)dt). Note that £ describes the strain- 
ing motions (e.g. Acheson 1990) which are divergence free and involve stretching and 
squashing in mutually perpendicular directions. 

There are three characteristic time scales that can be associated to a collapsing fluid 
element. The cooling time, t coo i (see appendix), and two dynamical time scales : T^ yn = 
—9/2 ± £ which correspond to the collapse time along the two principal axis of the 
fluid element. In the case of an isobaric evolution (see appendix for details), we have 
t coo i ~ Tdyn and the characteristic time of evolution of the fluid element will be the 
cooling time. If one takes the pressure gradient into account, the different time scales are 
not necessarily equal and more complex situations can arise. 

When there is no turbulence the kinetic energy of the incoming flow is mainly con- 
verted into thermal pressure at the shock while in the turbulent case it is converted both 
into thermal pressure and turbulent motion. Consequently, the shocked gas in the low 
turbulence simulation has a very short cooling time (because of their high pressure) and 
a very long dynamical time because most of the kinetic energy was turned into thermal 
pressure. On the contrary, in the highly turbulent simulation, the shocked gas has a much 
longer cooling time and a smaller dynamical time. 

In particular, if a dynamical time scale is shorter than the cooling time, it means that 
an axis is collapsing faster than the gas can cool. Therefore the pressure and the pressure 
gradient will rise along this particular direction which will eventually bounce and start 
to expand. 

In view of this, it is possible to identify two physical mechanisms which explain why 
thermally unstable gas is generated by turbulence. 

First, as illustrated by Fig. [3] and [71 the thermal pressure is higher in the case of a 
weakly turbulent flow than in the case of a very turbulent flow making the cooling time 
much shorter in the first case than in the second one. 
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The second mechanism which enhances the fraction of thermally unstable gas is a 
dynamical one related to the straining motions. Let us consider for reference an initially 
spheroidal piece of thermally unstable gas which contracts isotropically. In such case 9 is 
negative and E vanishes. The density, p, increases and because of the cooling, the internal 
energy decreases. Therefore if the dynamical time scale is larger than the cooling time, 
the pressure P decreases as well and the piece of gas keeps contracting. In this case the 
dynamical times are equal on both axis. Let us now consider the same piece of gas but 
with a positive E. The dynamical time scale will be reduced on the y-axis and enhanced 
on the x-axis. The cooling time will remain the same (at least in a first phase) since the 
evolution of the density and of the pressure depends only on 9. However, since the rate 
of contraction along the y-axis is higher, the pressure gradient will also be greater and 
the contraction can be slowed or even turned into expansion. 

All these dynamical effect are amplified by the thermally unstable nature of the gas. 
If for any dynamical reason the contraction is reduced, then the density will be lower 
and since d p P < the pressure will be higher and it will therefore be even more difficult 
to collapse. 

This stabilisation of unstable gas is illustrated in Fig. where the evolution of 
a fluid element in the pressure-density diagram is displayed. The bottom panel is for 
(9, E) = (0,0), whereas the upper panel is for (9, E) = (0, 0.8/r coo i). Since we want to 
illustrate the dynamical effect of E, we have kept the initial thermal pressure identical 
for the two models; but as was mentioned earlier, in the simulation the initial thermal 
pressure is lower on average when E is large. It is seen that in the first case, the gas 
contracts rapidly with almost no oscillations and the thermal pressure decreases rapidly. 
In the second case the fluid element oscillated rapidly and the evolution is more isobaric. 
The fluid particle remains in the unstable domain about twice the time it spends in the 
unstable domain when E = 0. 

Finally, we stress that both mechanisms presented to explain the presence of unstable 
gas are related to the presence of turbulence. It is expected that filamentary structures 
will be produced by such motions. Since straining motions (i.e. E) are produced in a tur- 
bulent flow, this mechanism explains how turbulence may produce filaments of thermally 
unstable gas. 

4.1.2. Predictions of the semi-analytical model 

We now present the quantitative predictions of Eqs. (|7|)- (|10|l for the time spent in the 
thermally unstable domain. We consider a fluid element of warm gas after it has been 
shocked and we follow its evolution until it reaches a temperature of 200 K. We explore 
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Fig. 11. Time spent by the fluid particle in the thermally unstable domain (defined as 
200 K < T < 5000 K) during its contraction as a function of E. The 5 curves represent 
respectively (from top to bottom) a mach number of 1.2, 1.4, 1.6, 1.8 and 2. 

a large ensemble of initial conditions defined by the initial state of the gas (temperature 
and density), the initial values of 9 and E as well as the size of the fluid particle, 

Fig. 1111 shows the mean time spent in the thermally unstable domain as a function 
of the initial value of E for 5 values of the initial mach number, M, namely 1.2, 1.4, 1.6, 
1.8 and 2. The initial density and temperature are obtained by applying the Rankinc- 
Hugoniot relations for an isothermal shock of mach number M and preshocked gas of 
density and temperature equal to the values used as initial conditions in our simulations 
i.e no = 0.76 cm -3 and To = 7100 K. For each value of M and S, we calculate the 
particle evolution for 6 ranging from -5.0 to 2.5 Myr -1 , for I ranging from 0.5 to 10 pc 
and compute the mean value. 

It is seen that, as expected, the time spent by the fluid particle in the thermally 
unstable domain, r„, decreases when the Mach number increases and increases when E 
increases. The ratio between these times for M = 2 and S = and M = 1.2 and E = 4 
is about 10. 

Under simple assumptions, the results presented in Fig. llll can be used to predict the 
relative ratio between the fraction of unstable and cold gas for different values of M and 
E. Let m u and m c be the mass of the unstable and cold gas respectively. The cold gas 
forms from the collapse of the unstable gas and disappears when it reaches the upper or 
lower sides of our computational box. Therefore one has: 

— - ~ — 11 

at t u t c 

where r c is the mean time spent by the cold gas in our simulation box. Therefore when 
statistical equilibrium is reached in our box, we have: 

t„(A/ 2 ,E 2 ) \m c J (MuSi) V m c/ (M2:E2) 
where we have assumed that r c does not vary significantly with M and E. 
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Fig. 12. Fraction of cold gas (full line) and fraction of thermally unstable gas plus cold 
gas (dashed line) as a function of < £ > for the 4 cases e =0.5, 2, 4 and 6 for various 
time steps. Each point correspond to the average over several outputs of a simulation 
with a given e and the error bars show the total dispersion of the results (100% of the 
points are within the error bar). 

4.2. Comparison with the numerical results 

In the model presented previously turbulence is able to stabilise transiently a piece of 
thermally unstable gas making it able to survive longer. This is done through the mecha- 
nisms presented in the previous section, which are both related to the presence of straining 
motions (i.e. £). Therefore this model predicts that the thermally unstable gas should 
be correlated with the parameter E in the numerical simulations. In order to verify this 
effect we have smoothed out the velocity field at a scale of 0.2 pc (10 pixels) and com- 
puted E at every pixel (Note that larger smoothing scales, e.g. 0.4 or 1 pc lead to similar 
results). Then the correlation between E and the fluid temperature has been investigated 
both globally and locally. 

Fig. E| displays the total fraction of cold gas (full line) and the total fraction of 
cold plus thermally unstable gas (dashed line) in the 4 simulations e = 0.5, 2, 4 and 
6 at different times as a function of < E >. It is seen that the second one is almost 
independent of < E > whereas the first decreases linearly with it. This result suggests 
that on one hand, the fraction of warm gas which is driven into the unstable regime 
depends mainly on the strength of the converging flow and is roughly independent of 
the level of turbulence. On the other hand, the thermally unstable gas is able to live 
longer when E is stronger. For e = 0.5, one has < E >~ 1.5 10~ 6 year -1 , m u /m c ~ 0.32 
whereas the mean pressure of the warm gas before it starts contracting in CNM is about 
3 10 -12 erg/cm 3 (obtained from Fig.[3J) corresponding to a mach number M, of ~ 1.8. 
For e = 4, < E >~ 3.5 10~ 6 , M ~ 1.2 (obtained from Fig. 01 and m u /m c ~ 2.8. 
The relation between the pressure and the Mach number are obtained using isothermal 
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Fig. 13. Probability distribution function (pdf) of E for e = 0.5 (dashed line), e = 2 
(dotted line), e = 4 (dashed-dotted line) and e — 6 (full line). Upper panel is for the 
whole simulation, middle panel for cold gas only and lower panel for thermally unstable 
gas. 

Rankine-Hugoniot relations, which seems a reasonable assumption looking at Fig. Eland 
[7| The ratio between m M /m c for these 2 values is therefore about 8. 

Let us make a quantitative comparison with the prediction of the semi-analytical 
model. Fig. ^| predicts for the case of the simulation with e = 0.5 (M ~ 1.8 and 
< E >~ 1.5 10~ 6 year" 1 ), r„ ~ 1.5 Myr whereas for e = 4 (M ~ 1.2 and < E >~ 3.5 10" 6 
year -1 ) it predicts about 10 Myr. According to the arguments presented in Sect 14. 1.21 
this leads to an estimate for the ratio between m u /m c in these 2 cases of about 10 / 1.5~ 
7 which is in very reasonable agreement with the value (~ 8) measured in the numerical 
experiment. 

Fig. El displays the pdf of E for e = 0.5 (dotted line) and e — 4 (full line). The upper 
panel is for the whole simulation whereas the middle one is restricted to cold gas (T < 
200 K) and the lower to thermally unstable gas (200 K < T < 5000 K). It is obvious 
that in the cold phase, E is much lower than for the thermally unstable gas which is 
associated with high values of E. Note that for e = 0.5, the thermally unstable gas has 
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a lower E than for e = 4. This is consistent with the fact that this is turbulence which 
generates flows having large value of S. 

Both results are, qualitatively and quantitatively, in good agreement with the expla- 
nation proposed in the previous sections. Another important, though qualitative, fact is 
that the thermally unstable gas is organised in filaments which is also a natural outcome 
of the mechanism based on straining motions. 

5. Conclusion 

We perform 2-d numerical simulations with the aim of understanding the dynamics of 
the thermally bistable interstellar atomic hydrogen. The size of our 1000 2 simulation box 
is 20 pc, leading to a numerical resolution of 0.02 pc. Such resolution is needed in order 
to properly describe the cold HI structures. In order to mimic a large scale converging 
flow that triggers the transition of the warm phase into the cold one, warm gas is injected 
continuously from 2 opposite sides of the box. The gas can freely escape at the 2 other 
boundaries. Various amplitudes of turbulent fluctuations have been applied. 

When the flow is weakly turbulent, a layer of compressed WNM forms and quickly 
fragments into structures of cold gas. When the flow is very turbulent, it is less organised. 
Nevertheless the thermally bistable behaviour is not erased and indeed locally very similar 
to the classical picture of the 2-phase medium. In particular, the 2 phases are connected 
through sharp thermal fronts and are locally in rough pressure equilibrium. 

In order to characterise the CNM structures found in the numerical simulation, wc 
apply a simple threshold algorithm to identify them and give the pdf of some of their 
intrinsic properties as mean density, temperature, velocity dispersion, size and aspect 
ratio, r. Typical values for these parameters are respectively: < n >~ 50 cm~ 3 , T~80 
K, Sv ~ 0.5 km/s, 1~ 0.1 pc and r~0.5. 

A large fraction of thermally unstable gas which increases with the amplitude of 
the turbulent forcing, is found. This thermally unstable gas tends to be organised in 
filamentary structures. 

A semi-analytical model for the thermal and dynamical evolution of a fluid particle 
is presented. It is based on the Lagrangian description of a fluid element. This model 
predicts that substantial straining motions may efficiently stabilise a piece of thermally 
unstable gas allowing it to survive a longer period of time. We then verify than both 
locally and globally the thermally unstable gas in the numerical simulations is indeed 
very strongly correlated with substantial straining motions. 
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Appendix A: A semi-analytical model 

In this appendix we derive the semi- analytical model that we use in Sect.|2| The idea of 
the model is to follow a fluid element by writing Eas. lllHl in a Lagrangian form and then, 
making some approximations, to compute the pressure gradients. 

The conservation of the mass can be written in a Lagrangian form as: 

p + P 9 = Q. (A.l) 

where 9 is the divergence of the velocity field and the dot represent a total (i.e. 
Lagrangian) derivative. The evolution of thermal energy can also be simply written in a 
Lagrangian form: 

£ + -e = -C(p,T), (A.2) 
P P 

where E — P/(p{ r f — 1)) is the specific thermal energy. 

In order to have an evolution equation for 9 we take the gradient of Eq. J2J. As it 
is common in fluid mechanics, we then decompose the resulting tensorial equation into 
its trace, symmetric trace- free and anti-symmetric part. For the sake of simplicity and 
since we want to compare our model to the simulation presented above, we write these 
equations for a bi-dimcnsional flow : 



9+ 9 - + 2(af + <j\ - uj 2 ) = -{P xx + P yy ), (A.3) 

01+ 9UI = -~(P XX - Pyy), (A.4) 

0=2+ 0(T 2 =-\(Pxy + Pjx), (A.5) 

lo+ 9u + uj(a 2 - <ti) =0, (A. 6) 

where o\ = 0.5(d x u x — d y u y ) and oi = 0.5(d x u y + d y u x ) are the two components of 
the shear tensor and uj = 0.5(d x u y — d y u x ) is the rotational. The Py are defined by : 
P l3 =d l {d, J P/ P ). 

The system (|A.l|) - (IA.6ll is exact but is not closed since there is no evolution equa- 
tion for Pij . Therefore it cannot be integrated unless one makes some approximation to 
compute the Pij. 

The simplest case to consider, as a starting point, is the evolution of an isobaric 
fluctuation. In that case P^ — and the evolution of the fluid element is given by: 

9 = — = — C {p, T) = = (A.7) 

P jPo cti cr 2 
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where Pq is the constant pressure and eo = fo/(7 — 1) is the corresponding internal 
energy. The cooling time, t coo i = —eo/C(p,T) is the characteristic time of evolution of 
the system. 

In order to treat the pressure gradient let us first simplify the previous system. If 
rotation is neglected, the shear tensor can be diagonalized and the eigenvector basis will 
not rotate during the evolution. In this eigenframe, the only component of the shear 
tensor is given by S = \J o\ + a\ and the fluid element can be viewed as an ellipsoid 
defined by its semi-major and minor axis of length a and b. The equation of evolution 
for E is identical to that of <j\ (|A.4|I with the pressure gradient taken in the eigenframe. 
The evolution of the ellipsoid axis is given by: 



a(t) = a Q exp (J (6/2 + Y,)dt 



a a e az , b(t) = b a exp (J (9/2 - Y;)dt\ = 6 ae/as 

(A.8) 

where we have defined ag = exp(J * 0/2dt) and as — exp(J Q * Sdi). 

If we further assume that the pressure and density inside and outside the ellipsoid 
are uniform then we may write: 

p _ -Pcxt - P(t) P C xt ~ P(t) p ~ n (A a) 

xx ~ p(t)a(tf ' vv ~ p(t)b(tf a d xy ~ ^ { yj 

Using this approximation for , it is now possible to integrate the system <|A.1|I - <|A.6|) 
to determine the evolution of the fluid element. 
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